Practical 1B - microbiota data

2024 NUTRIM microbiome & metabolome workshop

1 Intro

Now we’re going to look at the microbiota data, from the same Papa et al. 2012 study on paediatric IBD patients and controls.

A note on sequence data processing

These data have already been processed into a table: counts of how often each sequence occurs in each sample.

It started as huge fastq text files, full of As, Cs, Ts and Gs, but we will not practice sequence data processing today, because it takes quite a long time to run.

If you generate your data with us at Medical Microbiology, we will run the latest processing approaches on our computational infrastructure and provide you with the processed data.

Nowadays we do amplicon sequencing with Illumina MiSeq, HiSeq, or similar technologies, and denoise the output with DADA2 to produce ASV count tables.

The example data we will use are a little older. The amplicons were sequenced with 454 pyrosequencing and clustered into OTUs, but the core principles of analysis remain the same.

And as you will see later today, the same approaches can also be applied to taxonomic abundance count tables obtained from shotgun metagenomic sequencing.

First we’re going to use the R skills we practised in part 1A to inspect this data.

SOLUTION: (click me?)

• In this practical, solution blocks like this are available…
• But try and write code independently before looking at the answer!
• Attempting to recall knowledge or solve problems is proven to enhance learning.
• So don’t look unless you’re really stuck! 💪

After that, we’ll take a look at some specialist R packages for microbiome analysis.

2 Learning goals 🧠

• Practice general R skills to inspect microbiota count data and taxonomy
• Gain basic familiarity with specialist microbiome R packages: phyloseq and microViz
• Learn to calculate, visualise, and interpret standard richness and diversity indices: Observed Richness and (Effective) Shannon diversity

3 Load R packages 📦

library(readxl)
library(here)
library(tidyverse)

4 Read and inspect data 🔍

Read the metadata file from part 1A. We will read the RDS file "all_metadata.rds"

meta <- read_rds(file = here("data/papa2012/processed/all_metadata.rds"))

Read the count table: this is stored as a TSV (tab-separated variables) formatted text file.

counts <- read_tsv(file = here("data/papa2012/papa2012_OTU_count_table.txt"))

Rows: 91 Columns: 36350
── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr     (1): sample
dbl (36349): OTU_00001, OTU_00002, OTU_00003, OTU_00004, OTU_00005, OTU_0000...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

counts

# A tibble: 91 × 36,350
   sample OTU_00001 OTU_00002 OTU_00003 OTU_00004 OTU_00005 OTU_00006 OTU_00007
   <chr>      <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
 1 099_AX         0        56        30         0         0       103         0
 2 199_AX        29         0         9        17         4         5        57
 3 062_BZ         0        95        33         0         0         0        32
 4 194_AZ        75        49       228         0        30        35         0
 5 166_AX       138         3        12         0        13        88       103
 6 219_AX         0       347         2         2         0         0         1
 7 132_AX       121         0         7         0         0         0       181
 8 026_AX         0         0         0         0         0         0         0
 9 102_AZ         0         0         0         0         0         0         0
10 140_AX        47         0        16         1         1        46         4
# ℹ 81 more rows
# ℹ 36,342 more variables: OTU_00008 <dbl>, OTU_00009 <dbl>, OTU_00010 <dbl>,
#   OTU_00011 <dbl>, OTU_00012 <dbl>, OTU_00013 <dbl>, OTU_00014 <dbl>,
#   OTU_00015 <dbl>, OTU_00016 <dbl>, OTU_00017 <dbl>, OTU_00018 <dbl>,
#   OTU_00019 <dbl>, OTU_00020 <dbl>, OTU_00021 <dbl>, OTU_00022 <dbl>,
#   OTU_00023 <dbl>, OTU_00024 <dbl>, OTU_00025 <dbl>, OTU_00026 <dbl>,
#   OTU_00027 <dbl>, OTU_00028 <dbl>, OTU_00029 <dbl>, OTU_00030 <dbl>, …

Your first challenge:

Read the taxonomy table stored in "data/papa2012/papa2012_taxonomy_table.txt"

Call the object taxonomy

SOLUTION:

taxonomy <- read_tsv(file = here("data/papa2012/papa2012_taxonomy_table.txt"))

Rows: 36349 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr (7): OTU, Kingdom, Phylum, Class, Order, Family, Genus

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

taxonomy

# A tibble: 36,349 × 7
   OTU       Kingdom  Phylum         Class               Order      Family Genus
   <chr>     <chr>    <chr>          <chr>               <chr>      <chr>  <chr>
 1 OTU_00001 Bacteria Firmicutes     Clostridia          Clostridi… Rumin… Faec…
 2 OTU_00002 Bacteria Proteobacteria Gammaproteobacteria Enterobac… Enter… Esch…
 3 OTU_00003 Bacteria Firmicutes     Clostridia          Clostridi… Lachn… Blau…
 4 OTU_00004 Bacteria Firmicutes     Clostridia          Clostridi… Lachn… Rose…
 5 OTU_00005 Bacteria Bacteroidetes  Bacteroidia         Bacteroid… Bacte… Bact…
 6 OTU_00006 Bacteria Firmicutes     Clostridia          Clostridi… Rumin… Faec…
 7 OTU_00007 Bacteria Firmicutes     Clostridia          Clostridi… Lachn… Rose…
 8 OTU_00008 Bacteria Firmicutes     Clostridia          Clostridi… Rumin… Faec…
 9 OTU_00009 Bacteria Proteobacteria Gammaproteobacteria Enterobac… Enter… Esch…
10 OTU_00010 Bacteria Firmicutes     Negativicutes       Selenomon… Veill… Dial…
# ℹ 36,339 more rows

4.1 Taxonomy table

Now practice inspecting the taxonomy table by completing the following tasks:

Task 1

Check how many distinct genera there are.

Tip: use unique() and length()

SOLUTION:

taxonomy$Genus %>% unique() %>% length()

[1] 142

Task 2

How many OTUs are there in each Phylum?

Tip: use count() or table()

SOLUTION:

taxonomy %>% count(Phylum)

# A tibble: 9 × 2
  Phylum                        n
  <chr>                     <int>
1 Actinobacteria              307
2 Bacteroidetes             13109
3 Cyanobacteria/Chloroplast     1
4 Firmicutes                17653
5 Fusobacteria                275
6 Proteobacteria             4806
7 Synergistetes                 6
8 Verrucomicrobia              59
9 <NA>                        133

taxonomy$Phylum %>% table(useNA = "ifany")

.
           Actinobacteria             Bacteroidetes Cyanobacteria/Chloroplast 
                      307                     13109                         1 
               Firmicutes              Fusobacteria            Proteobacteria 
                    17653                       275                      4806 
            Synergistetes           Verrucomicrobia                      <NA> 
                        6                        59                       133 

Task 3

What genera are in the phylum Actinobacteria?

Tip: use filter()

SOLUTION:

taxonomy %>% filter(Phylum == "Actinobacteria") %>% count(Genus, sort = TRUE)

# A tibble: 9 × 2
  Genus               n
  <chr>           <int>
1 Bifidobacterium   170
2 Collinsella       111
3 Eggerthella         8
4 Actinomyces         7
5 <NA>                5
6 Gordonibacter       3
7 Asaccharobacter     1
8 Atopobium           1
9 Rothia              1

Or for just their names:

taxonomy %>% filter(Phylum == "Actinobacteria") %>% pull(Genus) %>% unique()

[1] "Collinsella"     "Bifidobacterium" "Eggerthella"     "Atopobium"      
[5] "Actinomyces"     NA                "Gordonibacter"   "Asaccharobacter"
[9] "Rothia"         

4.2 OTU count table

First let’s plot a histogram of OTU number 1.

counts %>% ggplot(aes(OTU_00001)) + geom_histogram(binwidth = 5)

Looks like there are a lot of zeros!

table(OTU00001_has_0_counts = counts$OTU_00001 == 0, useNA = "ifany")

OTU00001_has_0_counts
FALSE  TRUE 
   59    32 

---

Attempt the following tasks to explore further!

Task 1

Use filter() to plot only the non-zero entries for OTU 1.

Try also using + scale_x_log10() to transform the plot axis scale.

SOLUTION:

counts %>%
  filter(OTU_00001 != 0) %>%
  ggplot(aes(OTU_00001)) +
  geom_histogram(bins = 30) +
  scale_x_log10()

Task 2

Use mutate() to create a temporary log-transformed variable log_OTU1 and plot its distribution.

Note that you can’t do log(0), so you will need to add 1 to all values first, i.e. log10(OTU_00001 + 1)

SOLUTION:

counts %>%
  mutate(log_OTU1 = log10(OTU_00001 + 1)) %>%
  ggplot(aes(log_OTU1)) +
  geom_histogram(bins = 30)

Task 3

Plot OTU 1 against OTU 2 as a scatter plot using geom_point.

Remember to add 1 before log10 transformation of both variables.

SOLUTION:

counts %>%
  mutate(log_OTU1 = log10(OTU_00001 + 1)) %>%
  mutate(log_OTU2 = log10(OTU_00002 + 1)) %>%
  ggplot(aes(x = log_OTU1, y = log_OTU2)) +
  geom_point()

Task 4

But what type of bacteria do these amplicon sequences belong to?

Look up the classifications of OTU 1 and OTU 2 in the taxonomy table.

SOLUTION:

taxonomy %>% filter(OTU %in% c("OTU_00001", "OTU_00002"))

# A tibble: 2 × 7
  OTU       Kingdom  Phylum         Class               Order       Family Genus
  <chr>     <chr>    <chr>          <chr>               <chr>       <chr>  <chr>
1 OTU_00001 Bacteria Firmicutes     Clostridia          Clostridia… Rumin… Faec…
2 OTU_00002 Bacteria Proteobacteria Gammaproteobacteria Enterobact… Enter… Esch…

taxonomy %>%
  filter(OTU %in% c("OTU_00001", "OTU_00002")) %>%
  select(OTU, Phylum, Family, Genus) 

# A tibble: 2 × 4
  OTU       Phylum         Family             Genus               
  <chr>     <chr>          <chr>              <chr>               
1 OTU_00001 Firmicutes     Ruminococcaceae    Faecalibacterium    
2 OTU_00002 Proteobacteria Enterobacteriaceae Escherichia/Shigella

Task 5

Now make histograms and look up the taxonomy for the next thousand OTUs…

SOLUTION:

Okay, that was a joke. 🤡

It is possible to make a thousand plots, because R is very good at repetitive tasks.

But, this would not be very useful, because you could not look at them all.

In the next section of this practical, we will explore smarter ways to analyse microbiota compositions.

---

5 Assemble a phyloseq 🛠️

So far we have not used any R packages specialised for microbiome data.

Let’s do so, because it will make our next tasks a lot easier!

library(phyloseq)
library(microViz)

We will start by combining our three dataframes into one phyloseq object

5.1 OTU counts + taxonomy

phyloseq uses row or column names to match OTUs across the count and taxonomy tables.

In addition, the taxonomy table must contain only taxonomic ranks in descending order, and the count table must be converted to a numeric matrix.

count_matrix <- counts %>% select(-sample) %>% as.matrix()
rownames(count_matrix) <- counts$sample

tax_matrix <- taxonomy %>% select(!OTU) %>% as.matrix()
rownames(tax_matrix) <- taxonomy$OTU

ps <- phyloseq(
  otu_table(count_matrix, taxa_are_rows = FALSE), 
  tax_table(tax_matrix)
)

ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 36349 taxa and 91 samples ]
tax_table()   Taxonomy Table:    [ 36349 taxa by 6 taxonomic ranks ]

5.2 Adding sample metadata

phyloseq uses row names to match samples (across the OTU count table and the sample metadata).

head(sample_names(ps))

[1] "099_AX" "199_AX" "062_BZ" "194_AZ" "166_AX" "219_AX"

Notice that in the phyloseq object the sample names have underscores, but they have hyphens in the sample metadata.

meta

# A tibble: 90 × 17
   ID    sample case_control diagnosis activity active   gender ethnicity
   <chr> <chr>  <chr>        <fct>     <fct>    <fct>    <chr>  <chr>    
 1 099A  099-AX Case         UC        severe   active   female White    
 2 199A  199-AX Control      Other     control  control  female Other    
 3 062B  062-BZ Case         CD        mild     active   male   White    
 4 194A  194-AZ Case         UC        mild     active   male   White    
 5 166A  166-AX Case         UC        severe   active   female Black    
 6 219A  219-AX Case         UC        inactive inactive female <NA>     
 7 132A  132-AX Case         CD        mild     active   female White    
 8 026A  026-AX Case         UC        mild     active   male   White    
 9 102A  102-AZ Control      Other     control  control  male   White    
10 140A  140-AX Control      Other     control  control  female White    
# ℹ 80 more rows
# ℹ 9 more variables: family_history <chr>, age_years <dbl>,
#   immunosuppression_level <fct>, medication <chr>, ibd_fhx <lgl>, abx <lgl>,
#   steroids <lgl>, imsp_other <lgl>, imsp_any <lgl>

We can fix that with the stringr package str_replace function.

meta$sample <- meta$sample %>% str_replace(pattern = "-", replacement = "_")
meta

# A tibble: 90 × 17
   ID    sample case_control diagnosis activity active   gender ethnicity
   <chr> <chr>  <chr>        <fct>     <fct>    <fct>    <chr>  <chr>    
 1 099A  099_AX Case         UC        severe   active   female White    
 2 199A  199_AX Control      Other     control  control  female Other    
 3 062B  062_BZ Case         CD        mild     active   male   White    
 4 194A  194_AZ Case         UC        mild     active   male   White    
 5 166A  166_AX Case         UC        severe   active   female Black    
 6 219A  219_AX Case         UC        inactive inactive female <NA>     
 7 132A  132_AX Case         CD        mild     active   female White    
 8 026A  026_AX Case         UC        mild     active   male   White    
 9 102A  102_AZ Control      Other     control  control  male   White    
10 140A  140_AX Control      Other     control  control  female White    
# ℹ 80 more rows
# ℹ 9 more variables: family_history <chr>, age_years <dbl>,
#   immunosuppression_level <fct>, medication <chr>, ibd_fhx <lgl>, abx <lgl>,
#   steroids <lgl>, imsp_other <lgl>, imsp_any <lgl>

Now make them row names, check they all match, and add the metadata to the phyloseq object.

meta_df <- as.data.frame(meta)
rownames(meta_df) <- meta$sample
all(rownames(meta_df) %in% sample_names(ps))

[1] TRUE

sample_data(ps) <- meta_df
ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 36349 taxa and 90 samples ]
sample_data() Sample Data:       [ 90 samples by 17 sample variables ]
tax_table()   Taxonomy Table:    [ 36349 taxa by 6 taxonomic ranks ]

6 microViz 🦠👁️

microViz provides tools for working with phyloseq objects.

Let’s take a look with some basic microViz functions.

samdat_tbl(ps) # retrieve sample data as a tibble

# A tibble: 90 × 18
   .sample_name ID    sample case_control diagnosis activity active   gender
   <chr>        <chr> <chr>  <chr>        <fct>     <fct>    <fct>    <chr> 
 1 099_AX       099A  099_AX Case         UC        severe   active   female
 2 199_AX       199A  199_AX Control      Other     control  control  female
 3 062_BZ       062B  062_BZ Case         CD        mild     active   male  
 4 194_AZ       194A  194_AZ Case         UC        mild     active   male  
 5 166_AX       166A  166_AX Case         UC        severe   active   female
 6 219_AX       219A  219_AX Case         UC        inactive inactive female
 7 132_AX       132A  132_AX Case         CD        mild     active   female
 8 026_AX       026A  026_AX Case         UC        mild     active   male  
 9 102_AZ       102A  102_AZ Control      Other     control  control  male  
10 140_AX       140A  140_AX Control      Other     control  control  female
# ℹ 80 more rows
# ℹ 10 more variables: ethnicity <chr>, family_history <chr>, age_years <dbl>,
#   immunosuppression_level <fct>, medication <chr>, ibd_fhx <lgl>, abx <lgl>,
#   steroids <lgl>, imsp_other <lgl>, imsp_any <lgl>

# get the OTU table, or part of it
otu_get(ps, taxa = 1:5, samples = c("132_AX", "166_AX", "102_AZ"))

OTU Table:          [5 taxa and 3 samples]
                     taxa are columns
       OTU_00001 OTU_00002 OTU_00003 OTU_00004 OTU_00005
132_AX       121         0         7         0         0
166_AX       138         3        12         0        13
102_AZ         0         0         0         0         0

# get the taxonomy table
tt_get(ps) %>% head(3)

Taxonomy Table:     [3 taxa by 6 taxonomic ranks]:
          Kingdom    Phylum           Class                 Order              
OTU_00001 "Bacteria" "Firmicutes"     "Clostridia"          "Clostridiales"    
OTU_00002 "Bacteria" "Proteobacteria" "Gammaproteobacteria" "Enterobacteriales"
OTU_00003 "Bacteria" "Firmicutes"     "Clostridia"          "Clostridiales"    
          Family               Genus                 
OTU_00001 "Ruminococcaceae"    "Faecalibacterium"    
OTU_00002 "Enterobacteriaceae" "Escherichia/Shigella"
OTU_00003 "Lachnospiraceae"    "Blautia"             

7 Build some bar charts 📊

Sequencing data are compositional. The total number of reads per sample is mostly arbitrary, and so the counts should be interpreted as relative abundances instead of absolute abundances.

A simple way to visualise compositional data is as percentages, proportions of a whole.

Stacked bar charts are great for this: each bar represents one sample, and we can show the abundance data as proportions, after aggregating the counts by taxonomy.

ps %>% tax_fix() %>% comp_barplot("Phylum", n_taxa = 4, label = NULL)

This plot is aggregated into phyla, which is easy to read, and provides basic information.

We see that most samples contain a mixture of Firmicutes and Bacteroidetes, but some appear to be dominated by Proteobacteria instead.

Beware Fickle Phyla!

Beware, phylum names have changed recently!

• Actinobacteria is now Actinomycetota
• Bacteroidetes is now Bacteroidota
• Proteobacteria is now Pseudomonadota (!)
• Firmicutes is now Bacillota (!!)

In all prior research, you will see the original names, but in coming years, the new names will be increasingly adopted.

The best source for checking/searching official and alternative names is probably LPSN at bacterio.net

ps %>% tax_fix() %>% comp_barplot("Genus", n_taxa = 11, merge_other = F, label = NULL)

This plot is aggregated into genera, which provides more detail, but is harder to read, and we cannot give every genus a distinct colour.

We see that many samples contain a relatively large proportion of Bacteroides or Faecalibacterium but there is quite some further variation!

tax_fix() ?

Filling gaps in the taxonomy table

• You might have noticed that the taxonomy table has some NA values.
• This often occurs when a sequence cannot be uniquely classified at the Genus level.
• The short 16S amplicon sequenced may only allow a unique classification at Family rank, or above.

tt_get(ps) %>% as.data.frame() %>% filter(is.na(Genus)) %>% head(3)

           Kingdom         Phylum               Class             Order
OTU_00014 Bacteria     Firmicutes          Clostridia     Clostridiales
OTU_00102 Bacteria Proteobacteria Gammaproteobacteria Enterobacteriales
OTU_00151 Bacteria           <NA>                <NA>              <NA>
                      Family Genus
OTU_00014    Lachnospiraceae  <NA>
OTU_00102 Enterobacteriaceae  <NA>
OTU_00151               <NA>  <NA>

• We need to fill those gaps! and we can do this with tax_fix(), which copies info down from a higher rank to fill the gaps.
• It often works fine with default settings, but sometimes there are more complicated taxonomy table problems
• To look at and fix the taxonomy table interactively, try tax_fix_interactive(ps) in the console
• (You may need to allow popups on your browser!)

Press Stop!

Running tax_fix_interactive(ps) will open a new web browser tab.

If you don’t see anything after running the command, you might need to unblock pop-ups!

When you are done looking, click the red STOP 🛑 button in the R console!

Let’s update our ps phyloseq object with the default fix.

ps <- tax_fix(ps)

Check the first few taxa now look correct?

#

8 Discover diversity 🌳🌲🌴

As a last task for this introductory session, we will calculate and visualise alpha diversity.

Why is diversity interesting?

Biologically

• Higher diversity ecosystems are probably more resilient to perturbations
• Lower gut microbiota diversity sometimes accompanies various health problems (in adults)
• BUT: diverse == healthy is not TRUE for all ecosystems (e.g. early infant gut microbiome)
• So, consider your own data and diversity hypotheses carefully

Practically

• Diversity indices provide a simple “one number” summary of each ecosystem
• This makes it relatively easy to compare samples, and do statistical testing

8.1 Observed richness

• The simplest richness measure is just counting, a.k.a. “Observed Richness”.
• Let’s compute the observed richness of genera.
• ps_calc_richness() computes the index for each sample and adds it to your sample_data

ps_alpha <- ps %>%
  tax_fix() %>%
  ps_calc_richness(rank = "Genus", index = "observed", varname = "N_genera")

8.2 Shannon diversity

Next we will calculate the Shannon diversity index at the level of Genus.

• The Shannon index is a commonly used diversity measure, with this formula: \(H = -\sum_{i=1}^Np_i\ln p_i\)
• Shannon index is often is denoted with \(H\), and here \(p_i\) denotes the proportional abundance of the \(i\)’th of \(N\) taxa in the sample.

Explanation of the Shannon index formula

• For each taxon \(i\), you multiply its proportional abundance \(p_i\) by the natural log of that proportion \(\ln p_i\), and sum those values.
• Try it out for yourself to convince yourself you get larger (negative) values for higher proportions.
• The highest value you can achieve with \(N\) taxa occurs with equal proportions (e.g. with 20 taxa, maximum diversity occurs if each has a relative abundance of 5%, i.e. 0.05)
• Lastly, you change the sign of the result to a positive number, for ease of interpretation (this just makes more intuitive sense: as higher positive numbers indicates higher diversity.)

We will also compute an exponentiated version, the effective number of genera, or effective Shannon index.

Explanation of the Effective Shannon Diversity

The exponent of the Shannon index \(e^H\) represents the number of taxa (genera) that would be present in an evenly abundant ecosystem with the same Shannon index.

• The numeric value of the Shannon index itself has no intuitive meaning.
• You can compare them, but can’t easily interpret any one number.
• So, the concept of “effective numbers” of taxa is useful here.
• If your original ecosystem was actually perfectly even, then \(e^H = N\)
• Where N is the observed richness
• The more uneven an ecosystem, the further \(e^H\) will be from \(N\)

ps_alpha <- ps_alpha %>% 
  ps_calc_diversity(index = "shannon", rank = "Genus", varname = "Shannon_Genus") %>% 
  ps_mutate(Effective_Shannon_Genus = exp(Shannon_Genus))

ps_alpha

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 36349 taxa and 90 samples ]
sample_data() Sample Data:       [ 90 samples by 20 sample variables ]
tax_table()   Taxonomy Table:    [ 36349 taxa by 6 taxonomic ranks ]

The new genus-level Shannon diversity variables are stored in the sample data slot of ps_alpha

ps_alpha %>% samdat_tbl() %>% select(sample, Shannon_Genus, Effective_Shannon_Genus)

# A tibble: 90 × 3
   sample Shannon_Genus Effective_Shannon_Genus
   <chr>          <dbl>                   <dbl>
 1 099_AX          1.89                    6.63
 2 199_AX          2.71                   15.1 
 3 062_BZ          1.66                    5.25
 4 194_AZ          1.52                    4.59
 5 166_AX          1.92                    6.81
 6 219_AX          1.22                    3.39
 7 132_AX          2.35                   10.5 
 8 026_AX          2.14                    8.47
 9 102_AZ          2.56                   12.9 
10 140_AX          1.32                    3.74
# ℹ 80 more rows

8.3 Diversity Distributions

First we will plot simple histograms of the richness and diversity values we observe.
Notice the different range of values for each one.

Observed Richness

ps_alpha %>% samdat_tbl() %>% pull(N_genera) %>% hist(main = "Observed Richness, Genus")

Shannon \(H\)

ps_alpha %>% samdat_tbl() %>% pull(Shannon_Genus) %>% hist(main = "Shannon Diversity, Genus")

Effective Shannon \(e^H\)

ps_alpha %>% samdat_tbl() %>% pull(Effective_Shannon_Genus) %>% hist(main = "Effective Number of Genera")

8.4 Interpreting diversity

Observed Richness: Number of Genera

Each sample is sorted and labelled by the observed richness of genera

Can you spot samples with equal richness but clear differences in evenness?

ps_alpha %>% 
  ps_arrange(N_genera) %>% 
  comp_barplot(
    tax_level = "Genus", n_taxa = 19, merge_other = FALSE, 
    sample_order = "asis", label = "N_genera"
  ) +
  coord_flip()

Don’t forget to look at the other tab, where the bars are labelled by their diversity!

Diversity: Shannon Effective Number of Genera

Now each sample is labelled with the effective Shannon diversity of genera (\(e^H\))

Do you see the general relationship of \(e^H\) with sample composition?

ps_alpha %>% 
  ps_arrange(Effective_Shannon_Genus) %>% 
  ps_mutate(short_shannon = formatC(Effective_Shannon_Genus, digits = 1, format = "f")) %>% 
  comp_barplot(
    tax_level = "Genus", n_taxa = 19, merge_other = FALSE, 
    sample_order = "asis", label = "short_shannon"
  ) +
  coord_flip()

8.5 Save our phyloseq

We have assembled a phyloseq object and calculated richness and diversity measures.

Let’s store the result of this processing, by writing the phyloseq object to a file.

We can use the saveRDS() function to do this.

saveRDS(ps_alpha, file = here("data", "papa2012", "processed", "papa12_phyloseq.rds"))

9 Next! ⏩

• Its time to stop exploring and take a break.
• In the next lecture you will learn more about the main approaches for microbiota data analysis.
• In the next practical session, we will come back to these data with a plan, an analysis plan!
• Click here: https://david-barnett.github.io/2024-NUTRIM-microbiome/practical-2/web/practical2A-instructions.html

10 Session info

sessioninfo::session_info()

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.0 (2024-04-24)
 os       macOS Sonoma 14.5
 system   aarch64, darwin20
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       Europe/Amsterdam
 date     2024-06-29
 pandoc   3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package          * version  date (UTC) lib source
 ade4               1.7-22   2023-02-06 [1] CRAN (R 4.4.0)
 ape                5.8      2024-04-11 [1] CRAN (R 4.4.0)
 Biobase            2.64.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 BiocGenerics       0.50.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 biomformat         1.32.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 Biostrings         2.72.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 bit                4.0.5    2022-11-15 [1] CRAN (R 4.4.0)
 bit64              4.0.5    2020-08-30 [1] CRAN (R 4.4.0)
 ca                 0.71.1   2020-01-24 [1] CRAN (R 4.4.0)
 cellranger         1.1.0    2016-07-27 [1] CRAN (R 4.4.0)
 cli                3.6.2    2023-12-11 [1] CRAN (R 4.4.0)
 cluster            2.1.6    2023-12-01 [1] CRAN (R 4.4.0)
 codetools          0.2-20   2024-03-31 [1] CRAN (R 4.4.0)
 colorspace         2.1-0    2023-01-23 [1] CRAN (R 4.4.0)
 crayon             1.5.2    2022-09-29 [1] CRAN (R 4.4.0)
 data.table         1.15.4   2024-03-30 [1] CRAN (R 4.4.0)
 digest             0.6.35   2024-03-11 [1] CRAN (R 4.4.0)
 dplyr            * 1.1.4    2023-11-17 [1] CRAN (R 4.4.0)
 evaluate           0.23     2023-11-01 [1] CRAN (R 4.4.0)
 fansi              1.0.6    2023-12-08 [1] CRAN (R 4.4.0)
 farver             2.1.2    2024-05-13 [1] CRAN (R 4.4.0)
 fastmap            1.2.0    2024-05-15 [1] CRAN (R 4.4.0)
 forcats          * 1.0.0    2023-01-29 [1] CRAN (R 4.4.0)
 foreach            1.5.2    2022-02-02 [1] CRAN (R 4.4.0)
 generics           0.1.3    2022-07-05 [1] CRAN (R 4.4.0)
 GenomeInfoDb       1.40.1   2024-05-24 [1] Bioconductor 3.19 (R 4.4.0)
 GenomeInfoDbData   1.2.12   2024-05-26 [1] Bioconductor
 ggplot2          * 3.5.1    2024-04-23 [1] CRAN (R 4.4.0)
 glue               1.7.0    2024-01-09 [1] CRAN (R 4.4.0)
 gtable             0.3.5    2024-04-22 [1] CRAN (R 4.4.0)
 here             * 1.0.1    2020-12-13 [1] CRAN (R 4.4.0)
 hms                1.1.3    2023-03-21 [1] CRAN (R 4.4.0)
 htmltools          0.5.8.1  2024-04-04 [1] CRAN (R 4.4.0)
 htmlwidgets        1.6.4    2023-12-06 [1] CRAN (R 4.4.0)
 httr               1.4.7    2023-08-15 [1] CRAN (R 4.4.0)
 igraph             2.0.3    2024-03-13 [1] CRAN (R 4.4.0)
 IRanges            2.38.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 iterators          1.0.14   2022-02-05 [1] CRAN (R 4.4.0)
 jsonlite           1.8.8    2023-12-04 [1] CRAN (R 4.4.0)
 knitr              1.47     2024-05-29 [1] CRAN (R 4.4.0)
 labeling           0.4.3    2023-08-29 [1] CRAN (R 4.4.0)
 lattice            0.22-6   2024-03-20 [1] CRAN (R 4.4.0)
 lifecycle          1.0.4    2023-11-07 [1] CRAN (R 4.4.0)
 lubridate        * 1.9.3    2023-09-27 [1] CRAN (R 4.4.0)
 magrittr           2.0.3    2022-03-30 [1] CRAN (R 4.4.0)
 MASS               7.3-60.2 2024-04-24 [1] local
 Matrix             1.7-0    2024-03-22 [1] CRAN (R 4.4.0)
 mgcv               1.9-1    2023-12-21 [1] CRAN (R 4.4.0)
 microbiome         1.26.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 microViz         * 0.12.3   2024-06-06 [1] local
 multtest           2.60.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 munsell            0.5.1    2024-04-01 [1] CRAN (R 4.4.0)
 nlme               3.1-164  2023-11-27 [1] CRAN (R 4.4.0)
 permute            0.9-7    2022-01-27 [1] CRAN (R 4.4.0)
 phyloseq         * 1.48.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 pillar             1.9.0    2023-03-22 [1] CRAN (R 4.4.0)
 pkgconfig          2.0.3    2019-09-22 [1] CRAN (R 4.4.0)
 plyr               1.8.9    2023-10-02 [1] CRAN (R 4.4.0)
 purrr            * 1.0.2    2023-08-10 [1] CRAN (R 4.4.0)
 R6                 2.5.1    2021-08-19 [1] CRAN (R 4.4.0)
 Rcpp               1.0.12   2024-01-09 [1] CRAN (R 4.4.0)
 readr            * 2.1.5    2024-01-10 [1] CRAN (R 4.4.0)
 readxl           * 1.4.3    2023-07-06 [1] CRAN (R 4.4.0)
 registry           0.5-1    2019-03-05 [1] CRAN (R 4.4.0)
 reshape2           1.4.4    2020-04-09 [1] CRAN (R 4.4.0)
 rhdf5              2.48.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 rhdf5filters       1.16.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 Rhdf5lib           1.26.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 rlang              1.1.4    2024-06-04 [1] CRAN (R 4.4.0)
 rmarkdown          2.27     2024-05-17 [1] CRAN (R 4.4.0)
 rprojroot          2.0.4    2023-11-05 [1] CRAN (R 4.4.0)
 rstudioapi         0.16.0   2024-03-24 [1] CRAN (R 4.4.0)
 Rtsne              0.17     2023-12-07 [1] CRAN (R 4.4.0)
 S4Vectors          0.42.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 scales             1.3.0    2023-11-28 [1] CRAN (R 4.4.0)
 seriation          1.5.5    2024-04-17 [1] CRAN (R 4.4.0)
 sessioninfo        1.2.2    2021-12-06 [1] CRAN (R 4.4.0)
 stringi            1.8.4    2024-05-06 [1] CRAN (R 4.4.0)
 stringr          * 1.5.1    2023-11-14 [1] CRAN (R 4.4.0)
 survival           3.5-8    2024-02-14 [1] CRAN (R 4.4.0)
 tibble           * 3.2.1    2023-03-20 [1] CRAN (R 4.4.0)
 tidyr            * 1.3.1    2024-01-24 [1] CRAN (R 4.4.0)
 tidyselect         1.2.1    2024-03-11 [1] CRAN (R 4.4.0)
 tidyverse        * 2.0.0    2023-02-22 [1] CRAN (R 4.4.0)
 timechange         0.3.0    2024-01-18 [1] CRAN (R 4.4.0)
 TSP                1.2-4    2023-04-04 [1] CRAN (R 4.4.0)
 tzdb               0.4.0    2023-05-12 [1] CRAN (R 4.4.0)
 UCSC.utils         1.0.0    2024-05-06 [1] Bioconductor 3.19 (R 4.4.0)
 utf8               1.2.4    2023-10-22 [1] CRAN (R 4.4.0)
 vctrs              0.6.5    2023-12-01 [1] CRAN (R 4.4.0)
 vegan              2.6-6.1  2024-05-21 [1] CRAN (R 4.4.0)
 vroom              1.6.5    2023-12-05 [1] CRAN (R 4.4.0)
 withr              3.0.0    2024-01-16 [1] CRAN (R 4.4.0)
 xfun               0.44     2024-05-15 [1] CRAN (R 4.4.0)
 XVector            0.44.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)
 yaml               2.3.8    2023-12-11 [1] CRAN (R 4.4.0)
 zlibbioc           1.50.0   2024-04-30 [1] Bioconductor 3.19 (R 4.4.0)

 [1] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library

──────────────────────────────────────────────────────────────────────────────